SUBROUTINE RDCOOR(NCID,LOPLEFT)
USE PARKIND1  ,ONLY : JPIM     ,JPRB,    JPRD
USE YOMHOOK   ,ONLY : LHOOK    ,DR_HOOK, JPHOOK
USE MPL_MODULE
#ifdef DOC
! (C) Copyright 2000- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *RDCOOR* - Reading netCDF file containing point coordinates

!     Purpose.
!     --------
!        initialization surface coordinates
!        writing coordinates to output files
!        reading the catchment mask

!**   Interface.
!     ----------
!        *CALL* *RDCOOR*

!     Explicit arguments :
!     --------------------
!     INPUT
!        NCID      integer         NETCDF FILE ID
!         LOPLEFT: FLAG INDICATING WHETHER MORE POINTS NEED TO BE DONE

!     Implicit arguments :
!     --------------------


!     Method.
!     -------
!       Opens a file called 'surfclim' to read relevant fields


!     Externals.
!     ----------
!       NETCDF-utilities

!     Reference.
!     ----------

!     Author.
!     -------
!        Bart vd Hurk, KNMI

!     Modifications.
!     --------------
!        Original : 2000-07-13

!     ------------------------------------------------------------------

#endif
USE YOMDPHY  , ONLY : NCSS,NLON, NLAT, NLALO, NPOI,NVHILO,NTILES &
           &,DZLAT,DZLON,NCOOR,NPDONE
USE YOMCST   , ONLY : RPI
USE YOMGC1S  , ONLY : GELAM  ,GELAT     ,LMASK
USE YOMLUN1S , ONLY : NPOSGG ,NPOSEFL   ,NPOSRES   ,NPOSWAT &
                     &,NPOSSUS,NPOSSUB   ,NPOSEVA   ,NPOSCLD &
                     &,NULOUT, NPOSCLM &
                     &,NPOSCO2 ,NPOSVEG  ,NPOSEXT  ,NPOSBIO &
                     &,NPOSTIL  ,NPOSVTY   ,NPOSRES 

USE YOMLOG1S , ONLY : NDIMCDF &
           &,LWREFL   ,LWRWAT   ,LWRSUB   ,LWRSUS   ,LWREVA &
           &,LWRCLD   ,LWRGG,  LWRCLM &
           &,LWRCO2   ,LWRVEG   ,LWREXT &
           &,LWRBIO   ,LWRTIL   ,LWRVTY &
           &,NDIMCDF

USE YOESOIL1S, ONLY : RDAW
USE NETCDF
USE NETCDF_UTILS, ONLY: NCERROR
IMPLICIT NONE

!* Arguments
LOGICAL,INTENT(OUT) :: LOPLEFT     ! FLAG FOR RUNNING MORE GRID POINTS
INTEGER(KIND=JPIM),INTENT(IN) :: NCID

!* Netcdf interface
INTEGER NIDLON,NILON,NIDLAT,NILAT,IERR,NVATTS,NDIMS,NPOS,NVARID,&
     &   NVARTYP,NVDIMS,NVARS,NGATTS,IRECDIM
INTEGER ISTART1(1),ICOUNT1(1), ISTART2(2),ICOUNT2(2)
CHARACTER*100 CNAME
INTEGER,ALLOCATABLE :: KMASKX(:)
INTEGER,ALLOCATABLE :: KMASKY(:)
REAL*4,ALLOCATABLE :: ZLAT(:)
REAL*4,ALLOCATABLE :: ZLON(:)
REAL*4,ALLOCATABLE :: ZREAL(:)
REAL*8,ALLOCATABLE :: DZREAL(:)
REAL*8,ALLOCATABLE :: ZBOT(:)

INTEGER,ALLOCATABLE :: IVT(:),ITILE(:)

!* Local variables
INTEGER(KIND=JPRB), PARAMETER :: JPNCDF=15   !9 + 6(CTESSEL)
REAL(KIND=JPRB),ALLOCATABLE :: ZGELAT(:)
REAL(KIND=JPRB),ALLOCATABLE :: ZGELAM(:)
INTEGER(KIND=JPIM) :: IPOS(JPNCDF)
LOGICAL LPOS(JPNCDF),LOPRES
LOGICAL,ALLOCATABLE :: LMASKG(:)
INTEGER(KIND=JPIM) :: NX,NY,NPOIGLB,IX1,IX2,IY1,IY2,IX,IY,JL,IST,IEN,&
     &                ICT,JLA,JLO,NMX,NMY,J,JK,JVT,JTILE
REAL(KIND=JPRB) :: ZZZ
INTEGER(KIND=JPIM) :: IVAR,JD,DIMLEN,VARID
INTEGER(KIND=JPIM) :: MYPROC, NPROC
REAL(KIND=JPHOOK)    :: ZHOOK_HANDLE

#include "minmax.intfb.h"

IF (LHOOK) CALL DR_HOOK('RDCOOR',0,ZHOOK_HANDLE)

MYPROC = MPL_MYRANK()
NPROC  = MPL_NPROC()

!* check the dimensions
NILON=-999
NILAT=-999
IF( MYPROC == 1 ) THEN
  CALL NCERROR (NF90_INQUIRE(NCID, nDimensions=NDIMS),'FINDING N DIMS')
  DO JD=1,NDIMS
    CALL NCERROR(NF90_INQUIRE_DIMENSION(NCID, JD, name=CNAME, len=DIMLEN),'INQ DIMENSION')
    SELECT CASE (TRIM(CNAME))
      CASE('x','lon')
        NILON=DIMLEN
        WRITE(*,*)'DIMENSION x',TRIM(CNAME),' FOUND WITH LEN=',NILON
     CASE('y','lat')
        NILAT=DIMLEN
        WRITE(*,*)'DIMENSION y',TRIM(CNAME),' FOUND WITH LEN=',NILAT
    END SELECT
  END DO
ENDIF

CALL MPL_BROADCAST(NDIMS,KROOT=1,KTAG=100,CDSTRING='NDIMS')
CALL MPL_BROADCAST(NILON,KROOT=1,KTAG=200,CDSTRING='NILON')
CALL MPL_BROADCAST(NILAT,KROOT=1,KTAG=300,CDSTRING='NILAT')

IF (NDIMCDF == 1 ) THEN
  NILAT = 0
ENDIF

IF(NILON /= NLON .OR. NILAT /= NLAT) THEN
  WRITE(NULOUT,*)'NLON OR NLAT NOT SPECIFIED CORRECTLY:'
  WRITE(NULOUT,*)'NLON IN NAMELIST AND SURFCLIM: ',NLON,NILON
  WRITE(NULOUT,*)'NLAT IN NAMELIST AND SURFCLIM: ',NLAT,NILAT
  CALL ABORT
ENDIF

!* read coordinates

IF(NDIMCDF == 2)THEN
  NX=NLON
  NY=NLAT
ELSE
  NX=NLALO
  NY=NLALO
ENDIF

ALLOCATE (ZLON(NX))
ALLOCATE (ZLAT(NY))
ALLOCATE (ZREAL(NLALO))
ALLOCATE (DZREAL(NLALO))
ALLOCATE (LMASKG(NLALO))
ALLOCATE (ZGELAT(NLALO))
ALLOCATE (ZGELAM(NLALO))
ALLOCATE (DZLON(NX))
ALLOCATE (DZLAT(NY))
ALLOCATE (LMASK(NLALO))
ALLOCATE (ZBOT(NCSS))

ZBOT(1)=0.5*RDAW(1)
DO JK=2,NCSS
  ZBOT(JK)=ZBOT(JK-1)+0.5*(RDAW(JK-1)+RDAW(JK))
ENDDO

ALLOCATE (IVT(NVHILO))
DO JVT=1,NVHILO
  IVT(JVT)=JVT
ENDDO
ALLOCATE (ITILE(NTILES))
DO JTILE=1,NTILES
  ITILE(JTILE)=JTILE
ENDDO  


IF( MYPROC == 1 ) THEN
  CALL NCERROR( NF90_INQ_VARID(NCID, 'lat', VARID),'getting varid lat' )
  CALL NCERROR( NF90_GET_VAR(NCID,VARID,DZLAT,(/1/),(/NY/)),'READING LAT')

  CALL NCERROR( NF90_INQ_VARID(NCID, 'lon', VARID),'getting varid lon' )
  CALL NCERROR( NF90_GET_VAR(NCID,VARID,DZLON,(/1/),(/NX/)),'READING LON')
ENDIF

CALL MPL_BROADCAST(DZLAT(:),KROOT=1,KTAG=100,CDSTRING='DZLAT')
CALL MPL_BROADCAST(DZLON(:),KROOT=1,KTAG=200,CDSTRING='DZLON')
!* -- catchment mask

IF(NDIMCDF == 2)THEN
  ISTART2(1)=1
  ISTART2(2)=1
  ICOUNT2(1)=NLON
  ICOUNT2(2)=NLAT
ELSE
  ISTART2(1)=1
  ICOUNT2(1)=NLALO
ENDIF

IF( MYPROC == 1 ) THEN
  CALL NCERROR( NF90_INQ_VARID(NCID, 'Mask', VARID),'getting varid Mask' )
  CALL NCERROR( NF90_GET_VAR(NCID,VARID,DZREAL,ISTART2,ICOUNT2),'READING Mask')
  LMASKG(:)=(ABS(DZREAL(:)-1.) < 0.0001)
ENDIF

CALL MPL_BROADCAST(LMASKG(:),KROOT=1,KTAG=100,CDSTRING='LMASKG')

!1s debugmode 
!LMASKG(1:NLALO)=.FALSE.
!LMASKG(44448)=.TRUE.   !tropical forest (central africa)
!LMASKG(9187)=.TRUE.   !berms_oa
!LMASKG(60500)=.TRUE.
IF(NCOOR.EQ.0)THEN
  NPOIGLB=COUNT(LMASKG)
ELSE
  NPOIGLB=1
ENDIF
LMASK=.FALSE.
!* extraction of points in sub-slabs:
!  also the coordinates of the slab must be stored for later writing
!  when 1D, NSLAB = nr of points                ! NSLAB DEPRECATED
!  when 2D, NSLAB = nr of latitude lines        ! NSLAB DEPRECATED
IF(NDIMCDF == 2)THEN
  IX1=1
  IX2=NLON
  IY1=1
  IY2=NLAT
  DO IX=IX1,IX2
    DO IY=IY1,IY2
      JL=IX+(IY-1)*NLON
      IF(NCOOR == 0.OR. JL == NCOOR)THEN
        LMASK(JL)=LMASKG(JL)
      ENDIF
    ENDDO
  ENDDO
  NPOI=COUNT(LMASK)
  IF(IY2 == NLAT)LOPLEFT=.FALSE.
ELSE
  IST=1
  IEN=NLALO
  DO JL=IST,IEN
    IF(NCOOR == 0.OR. JL == NCOOR)THEN
      LMASK(JL)=LMASKG(JL)
    ENDIF
  ENDDO
  NPOI=COUNT(LMASK)
  NPDONE=NPDONE+NPOI
  IF(NPDONE.GE.NPOIGLB)LOPLEFT=.FALSE.
ENDIF

!* nr of calculation grid points
WRITE(NULOUT,*)'TOTAL NUMBER OF GRID POINTS:            ',NLALO
WRITE(NULOUT,*)'NUMBER OF GRID POINTS IN CATCHMENT:     ',NPOIGLB
WRITE(NULOUT,*)'NUMBER OF GRID POINTS IN CURRENT BATCH: ',NPOI

!IF(NDIMCDF == 1)THEN
!* in case of 1D-grid, get the mapping coordinates
!  ALLOCATE (KMASKX(NLALO))
!  ALLOCATE (KMASKY(NLALO))
! ! check presence
!   LOPRES=.FALSE.
!   DO IVAR = 1,NVARS
!     CALL NCVINQ (NCID, IVAR, CNAME, NVARTYP,NVDIMS, NDIMS,NVATTS, IERR)
!     IF(CNAME(1:11).EQ.'Grid_Mask_x')THEN
!       LOPRES=.TRUE.
!       GOTO 101
!     ENDIF
!   ENDDO
! 101 CONTINUE
!   IF(LOPRES)THEN
!     NVARID = NCVID(NCID, 'Grid_Mask_x', IERR)
!     CALL NCVGT(NCID, NVARID, 1, NLALO, KMASKX, IERR)
!     NVARID = NCVID(NCID, 'Grid_Mask_y', IERR)
!     CALL NCVGT(NCID, NVARID, 1, NLALO, KMASKY, IERR)
!   ELSE
!    WRITE(NULOUT,*)'Grid_Mask_x not found'
!    KMASKX=0
!    KMASKY=0
!   ENDIF
!ENDIF

!* -- convert to radians

ZZZ=RPI/180.
IF(NDIMCDF == 2)THEN
  ICT=0
  DO JLA=1,NLAT
    DO JLO=1,NLON
      ICT=ICT+1
      ZGELAT(ICT)=DZLAT(JLA)*ZZZ
      ZGELAM(ICT)=DZLON(JLO)*ZZZ
      IF (ZGELAM(ICT) < 0.) THEN
        ZGELAM(ICT)=ZGELAM(ICT)+2.*RPI
      ENDIF
    ENDDO
  ENDDO
  NMX=NLON
  NMY=NLAT
ELSE
  DO JL=1,NLALO
    ZGELAT(JL)=DZLAT(JL)*ZZZ
    ZGELAM(JL)=DZLON(JL)*ZZZ
    IF (ZGELAM(JL) < 0.) THEN
      ZGELAM(JL)=ZGELAM(JL)+2.*RPI
    ENDIF
  ENDDO
  NMX=NLALO
  NMY=1
ENDIF

IF( MYPROC == 1 ) THEN
  CALL MINMAX('LATITUDE',ZGELAT,NMX,NMY,LMASK,NULOUT)
  CALL MINMAX('LONGITUD',ZGELAM,NMX,NMY,LMASK,NULOUT)
ENDIF

!* -- Pack the arrays

ALLOCATE (GELAT(NPOI))
ALLOCATE (GELAM(NPOI))
GELAT=PACK(ZGELAT,LMASK)
GELAM=PACK(ZGELAM,LMASK)

DEALLOCATE (ZLON)
DEALLOCATE (ZLAT)
DEALLOCATE (ZREAL)
DEALLOCATE (DZREAL)
DEALLOCATE (LMASKG)
DEALLOCATE (ZGELAT)
DEALLOCATE (ZGELAM)


!IF(NDIMCDF == 1)THEN
!  DEALLOCATE (KMASKX)
!  DEALLOCATE (KMASKY)
!ENDIF
WRITE(NULOUT,*) 'RDCOOR is over'

CALL MPL_BARRIER()

IF (LHOOK) CALL DR_HOOK('RDCOOR',1,ZHOOK_HANDLE)

END SUBROUTINE RDCOOR
